To conduct an accessibility analysis, we need: 1) OSM Road Data 2) Shapefile of the study area 3) GTFS feed
We use these inputs to make queries on Open Trip Planner
Define the geographic extents as a vector of c(xmin, ymin, xmax, ymax)
library(osmdata)
library(tidyverse)
# opq used to build an overpass query
# we define the geographic extents of the area we want to query from OSM. I have added the extents of the GCR as a vector
# add_osm_feature() used to specify what we are querying. It takes key, value pairs but here I want all types of highways so I specify the key only
# For a list of specific values to query, see https://wiki.openstreetmap.org/wiki/Key:highway
query <- opq(bbox = c(30.801369, 29.705080, 31.874483, 30.390664)) %>%
add_osm_feature(key = 'highway')
# here we want the data in xml format. Ideally we would use pbf data, as it is more compact and makes the analysis faster, but osmdata_pdf() returns empty files. Issue discussed here https://github.com/ropensci/osmdata/issues/74
osmdata_xml(query, filename = 'greater-cairo.osm')
# the xml file is downloaded to the project working directory
You can plot to see if the roads have been downloaded. This takes around 10 minutes and the sp file is huge. I will not do it here, but for checking purposes, this is how it would be done:
convert xml to sp using osmdata_sp(query)
plot sp::plot(highways_greaterCairo$osm_lines) highways_greaterCairo <- osmdata_sp(query)
An alternative to xml is pbf files. These are more compact and presumably make the forthcoming analysis faster. They can be downloaded from the HOT OSM export tool https://export.hotosm.org/en/v3/exports/new/formats?fbclid=IwAR3m_1bDZK2sYsA-jtPRPYGg9R7kTqt5hzjP88x3p2yvRd4i9tr11i2tG60
I tried to use a pbf file but the analysis was not quicker so I stuck to the xml file extracted from r
In this step I import a shapefile of Greater Cairo. The shapefile has the region divided into hexagons with the population and number of jobs in each hexagon
library(sf)
Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
# import the shapefile
cairo_hexagons <- st_read("Cairo Shapefiles/H3-res8/H3res-8_GCR_4326.shp")
Reading layer `H3res-8_GCR_4326' from data source `/Users/Hussein/Documents/Code Repos/GIS-Coursework/Cairo Shapefiles/H3-res8/H3res-8_GCR_4326.shp' using driver `ESRI Shapefile'
Simple feature collection with 1613 features and 22 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 30.8684 ymin: 29.71863 xmax: 31.83956 ymax: 30.35436
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
It is imported and in the corrected crs EPSG:4326
Let’s plot to check if it looks right
library(tmap)
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
# plot
tm_shape(cairo_hexagons) +
tm_polygons()
It looks right but the area in the top right needs to be removed. This is the 10th of Ramadan City and is not part of the GCR
# Remove 10th of Ramadan (it is in Sharqeya)
# dplyr::filter with !grepl. This removes all rows that contain 10th of Ramadan OR '10 of Rammadan' in nam_tfc column
cairo_hexagons <- cairo_hexagons %>% filter(!grepl('10th of Ramadan|10 Of Rammadan', nam_tfc))
# plot to check that it worked
tm_shape(cairo_hexagons) +
tm_polygons()
The next step is to get the centroid of each hexagon. This is because Open Trip Planner takes queries from lat lon coordinates
# get centroids in seperate feature
h3_centroids <- st_centroid(cairo_hexagons)
st_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data
# centroids are provided as geometry but we need the lat and lon in two seperate columns
# I use this function to split c(lat, lon) to two seperate columns
# FROM JM London (https://github.com/r-spatial/sf/issues/231)
# lat = Y lon = X
sfc_as_cols <- function(x, names = c("lon","lat")) {
stopifnot(inherits(x,"sf") && inherits(sf::st_geometry(x),"sfc_POINT"))
ret <- sf::st_coordinates(x)
ret <- tibble::as_tibble(ret)
stopifnot(length(names) == ncol(ret))
x <- x[ , !names(x) %in% names]
ret <- setNames(ret,names)
dplyr::bind_cols(x,ret)
}
# add lon and lat columns to dataframe using sfc_as_cols function
h3_centroids <- sfc_as_cols(h3_centroids)
# two additional columns (lat and lon) have been added to the sfc
h3_centroids
Simple feature collection with 1507 features and 24 fields
geometry type: POINT
dimension: XY
bbox: xmin: 30.87314 ymin: 29.72468 xmax: 31.76349 ymax: 30.28386
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
First 10 features:
index city_nm prnt_idx W3W_id H3_arkm znng_t_
1 883f5b3601fffff 6th of October 873f5b360ffffff towels.dissolve.puppets 1175299 Outer
2 883f5b3603fffff 6th of October 873f5b360ffffff freshest.burden.enacted 1175322 Outer
3 883f5b3605fffff 6th of October 873f5b360ffffff repaid.glaze.dislodge 1175197 Outer
4 883f5b3607fffff 6th of October 873f5b360ffffff ocean.reason.mermaids 1175220 Outer
5 883f5b3609fffff 6th of October 873f5b360ffffff intrigues.garden.tent 1175379 Outer
6 883f5b360bfffff 6th of October 873f5b360ffffff builders.irritate.hockey 1175402 Outer
7 883f5b3611fffff 6th of October 873f5b361ffffff handover.mainly.bends 1175448 Outer
8 883f5b3613fffff 6th of October 873f5b361ffffff picturing.trip.underway 1175470 Outer
9 883f5b3615fffff 6th of October 873f5b361ffffff whirlwind.huddle.encoder 1175345 Outer
10 883f5b3617fffff 6th of October 873f5b361ffffff contained.hobbit.sandpaper 1175368 Outer
nam_tfc Uniq_ID SSc_N_1 n n_bank n_employ n_commer n_health n_recrea n_pub CBD Hus_ratio
1 Industrial Zone 3 89 <NA> 6 0 6 0 0 0 0 <NA> 0.6
2 Industrial Zone 3 89 <NA> 6 0 6 0 0 0 0 <NA> 0.6
3 Industrial Zone 3 Extenstion 550 <NA> 0 NA NA NA NA NA NA <NA> NA
4 Industrial Zone 3 Extenstion 550 <NA> 1 0 1 0 0 0 0 <NA> 0.1
5 Industrial Zone 3 89 <NA> 2 0 2 0 0 0 0 <NA> 0.1
6 Industrial Zone 2 285 <NA> 7 1 6 0 0 0 0 <NA> 0.6
7 Industrial Zone 1 544 <NA> 2 0 2 0 0 0 0 <NA> 0.2
8 Industrial Zone 4 475 <NA> 19 18 0 0 0 0 1 <NA> 7.8
9 Industrial Zone 2 285 <NA> 1 1 0 0 0 0 0 <NA> 0.1
10 Industrial Zone 2 285 <NA> 0 NA NA NA NA NA NA <NA> NA
pop2018cap jobsLFScou area popdenscap geometry lon lat
1 0 10281 1.175 0 POINT (30.90253 29.90795) 30.90253 29.90795
2 0 10281 1.175 0 POINT (30.89874 29.9175) 30.89874 29.91750
3 0 0 1.175 0 POINT (30.89731 29.89941) 30.89731 29.89941
4 0 1713 1.175 0 POINT (30.89353 29.90896) 30.89353 29.90896
5 0 1900 1.175 0 POINT (30.91154 29.90695) 30.91154 29.90695
6 0 10333 1.175 0 POINT (30.90775 29.91649) 30.90775 29.91649
7 0 5305 1.175 0 POINT (30.90017 29.93559) 30.90017 29.93559
8 0 1463 1.175 0 POINT (30.89638 29.94514) 30.89638 29.94514
9 0 52 1.175 0 POINT (30.89495 29.92705) 30.89495 29.92705
10 0 2524 1.175 0 POINT (30.89116 29.9366) 30.89116 29.93660
I use the package by Malcolm Morgan to query OTP through R https://docs.ropensci.org/opentripplanner/index.html
library(opentripplanner)
# OTP expects its data to be stored in a specific structure
# I create a new folder called Open-Trip-Planner in my wd
#dir.create() creates a subfolder called "OTP-Cairo-All"
path_data <- file.path("Open-Trip-Planner", "OTP-Cairo-All")
dir.create(path_data)
'Open-Trip-Planner/OTP-Cairo-All' already exists
Now we need to download OTP and save it to the “OTP-Cairo-All” subfolder
# otp_dl_jar function will download the OTP and save it in the folder we created
# The function returns the path to the OTP jar file.
path_otp <- otp_dl_jar(path_data)
The OTP will be saved to Open-Trip-Planner/OTP-Cairo-All/otp.jar
trying URL 'https://repo1.maven.org/maven2/org/opentripplanner/otp/1.4.0/otp-1.4.0-shaded.jar'
Content type 'application/java-archive' length 61838786 bytes (59.0 MB)
==================================================
downloaded 59.0 MB
The next step is to add the GTFS and OSM files to the created folder. Create a subfolder in OTP-Cairo-All called graphs then create a subfolder in graphs called default Add the OSM and GTFS inside ‘default’ I followed this folder structure:
I add a GTFS file with all public transport in the GCR (Formal + Informal)
This data is used to build a Graph object which is the base for OpenTripPlanner
# Building an OTP Graph
# This code will create a new file Graph.obj that will be saved in the location defined by path_data.
# memory argument assigns more memory to building graph (to speed it up). R assigns 2GB by default
log <- otp_build_graph(otp = path_otp, dir = path_data, memory = 6000, analyst = TRUE)
2020-01-09 09:13:05 Basic checks completed, building graph, this may take a few minutes
The graph will be saved to Open-Trip-Planner/OTP-Cairo-All
2020-01-09 09:15:13 Graph built
Now we are ready to launch OpenTripPlanner
Now that OTP is launched and connected to r, we can start querying
Lets plot one of the isochrones to see if this worked
tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(reachPolygons[[568]]) +
tm_fill(col = "antiquewhite2") +
tm_borders()
The shape reachPolygons[[568]] is invalid. See sf::st_is_validrow names were found from a short variable and have been discarded
Now let’s do the analysis again using a GTFS feed with formal transport modes only. This means we have to subset the GTFS feed and keep only the formal agencies. In OpenTripPlanner, there is a bannedAgencies argument, but this does not exist in the OpenTripPlanner R Package (see the documentation here: https://docs.ropensci.org/opentripplanner/reference/otp_isochrone.html?q=otp%20_%20isoch )
This is definitely a useful enhancement and I am thinking of requesting it on Github
There is one R package that handles GTFS feeds but it doesn’t filter by agencies https://github.com/ipeaGIT/gtfs2gps/blob/master/vignettes/intro_to_gtfs2gps.md
I ended up using a java application that can be accessed from the command line. It subsets GTFS feeds in a matter of seconds. Check the ‘How to Reduce your GTFS’ section http://developer.onebusaway.org/modules/onebusaway-gtfs-modules/1.3.4-SNAPSHOT/onebusaway-gtfs-transformer-cli.html
I added a text file with arguments for the agencies I want to retain (formal agencies) {“op”:“retain”, “match”:{“file”:“agency.txt”, “agency_id”:“CTA”}} {“op”:“retain”, “match”:{“file”:“agency.txt”, “agency_id”:“CTA_M”}} {“op”:“retain”, “match”:{“file”:“agency.txt”, “agency_id”:“NAT”}}
Now repeat the steps done above to create a directory
First thing is to disconnect from OTP. i found that if I don’t do this, OTP runs the analysis on the old graph that has the GTFS file with all agencies (even though I provide a new path…)
otp_stop(warn=FALSE)
Now let’s give it the new path
library(opentripplanner)
# create subfolder called "OTP-Cairo-Formal"
path_data <- file.path("Open-Trip-Planner", "OTP-Cairo-Formal")
dir.create(path_data)
'Open-Trip-Planner/OTP-Cairo-Formal' already exists
# otp_dl_jar function will download the OTP and save it in the folder we created
# The function returns the path to the OTP jar file.
path_otp <- otp_dl_jar(path_data)
The OTP will be saved to Open-Trip-Planner/OTP-Cairo-Formal/otp.jar
trying URL 'https://repo1.maven.org/maven2/org/opentripplanner/otp/1.4.0/otp-1.4.0-shaded.jar'
Content type 'application/java-archive' length 61838786 bytes (59.0 MB)
==================================================
downloaded 59.0 MB
Before building the graph, add the GTFS feed and the osm road layer in the correct folder. See 3.1.1: Setting UP OTP
Now we can run the analysis again using formal agencies only.
Now we’re done with otp. Make sure to terminate the Java OTP instance (it takes up 2GB of precious memory)
# warn equals false so that it doesn't prompt me to confirm
otp_stop(warn=FALSE)
The following Java instances have been found:
47374 ?? 5:02.56 /usr/bin/java -Xmx2048M -jar Open-Trip-Planner/OTP-Cairo-Formal/otp.jar --router default --graphs Open-Trip-Planner/OTP-Cairo-Formal/graphs --server --port 8080 --securePort 8081
47711 ?? 0:00.00 sh -c ps -A |grep java
47713 ?? 0:00.00 grep java
character(0)
Now that we have the isochrones, we need to calculate how many jobs each isochrone intersects with.
For each intersected hexagon, we get: (area of intersection / total area of hexagon) * jobs in Hexagon
We then sum all the results to get the number of jobs accessible for the hexagon we are querying from
ALL AGENCIES:
# empty list to store output
totalJobs <-list()
library(lwgeom) # for st_make_valid (this handles invalid geometries https://github.com/r-spatial/sf/issues/870)
for (i in 1:nrows){
totalJobs[[i]] <- 0 # there are some points that OTP couldn't route from. try() used to assign 0 value to these points
try({
totalJobs[[i]] <- reachPolygons[[i]] %>%
st_make_valid() %>% # some geometries are invlid (this prevents an error)
st_intersection(cairo_hexagons) %>% # intersect reachPolygon with cairo_hexagons
mutate(int_area = (st_area(.)/1000000) %>% as.numeric()) %>% # add column with intersection area with each hexagon
mutate(int_jobs = (int_area/area)*jobsLFScou) %>% # add column with int_jobs: jobs intersected
summarise(total_jobs = sum(int_jobs)) # summarize and get the sum of jobs intersected by reachPolygon
})
}
FORMAL AGENCIES:
# empty list to store output
totalJobsFormal <-list()
library(lwgeom) # for st_make_valid (this handles invalid geometries https://github.com/r-spatial/sf/issues/870)
for (i in 1:nrows){
totalJobsFormal[[i]] <- 0 # there are some points that OTP couldn't route from. try() used to assign 0 value to these points
try({
totalJobsFormal[[i]] <- reachPolygonsFormal[[i]] %>%
st_make_valid() %>% # some geometries are invlid (this prevents an error)
st_intersection(cairo_hexagons) %>% # intersect reachPolygon with cairo_hexagons
mutate(int_area = (st_area(.)/1000000) %>% as.numeric()) %>% # add column with intersection area with each hexagon
mutate(int_jobs = (int_area/area)*jobsLFScou) %>% # int_jobs is the jobs intersected
summarise(total_jobs = sum(int_jobs)) # gets one value for jobs intersected by reachPolygon
})
}
Now that we have the job reach of each hexagon, we need to transfer these values form the totalJobs list back to the cairo_hexagons sf so that we can calculate accessibility scores and, eventually, plot the results
# add a column for the number of jobs reachable within 60 minutes using ALL MODES of public transport
for (i in 1:nrows){
cairo_hexagons$jobs_60_all[i] <- 0 #set default value = O: will be given to bad geometries
try({
cairo_hexagons$jobs_60_all[i] = totalJobs[[i]][[1]] # add totalJobs to new column in cairo_hexagons called jobs_60_all
})
}
#add column for the number of jobs reachable within 60 minutes using FORMAL MODES public transport
for (i in 1:nrows){
cairo_hexagons$jobs_60_formal[i] <- 0 #set default value = O: will be given to bad geometries
try({
cairo_hexagons$jobs_60_formal[i] = totalJobsFormal[[i]][[1]] # add totalJobs to new column in cairo_hexagons called jobs_60_all
})
}
Now lets get the accessibility score for each hexagon. This is the % of the total jobs in the GCR that are reachable from this hexagon
# percentage of jobs accessible from each hexagon - ALL MODES OF PUBLIC TRANSPORT
cairo_hexagons$access_per_60_all = ((cairo_hexagons$jobs_60_all/ sum(cairo_hexagons$jobsLFScou))*100)
# percentage of jobs accessible from each hexagon - FORMAL MODES OF PUBLIC TRANSPORT
cairo_hexagons$access_per_60_formal = ((cairo_hexagons$jobs_60_formal/ sum(cairo_hexagons$jobsLFScou))*100)
To calculate the accessibility score for the entire study area:
1- weigh each hexagons job reach/accessibility score by its population (multiply them)
2- divide the sum by the total population of the GCR
# as.numeric used to return a number instead of a matrix
# Average jobs reached using all modes
GCR_avg_jobs_all <- as.numeric((cairo_hexagons$jobs_60_all %*% cairo_hexagons$pop2018cap) / sum(cairo_hexagons$pop2018cap))
cat(" Average Job Reach Using All Modes :", GCR_avg_jobs_all, "\n")
Average Job Reach Using All Modes : 276254.4
# Average accessibility using all Modes (%):
GCR_avg_acc_all <- as.numeric((cairo_hexagons$access_per_60_all %*% cairo_hexagons$pop2018cap) / sum(cairo_hexagons$pop2018cap))
cat(" Average Accessibility Using All Modes :", GCR_avg_acc_all, "\n")
Average Accessibility Using All Modes : 4.71886
# Average jobs reached using formal modes
GCR_avg_jobs_formal <- as.numeric((cairo_hexagons$jobs_60_formal %*% cairo_hexagons$pop2018cap) / sum(cairo_hexagons$pop2018cap))
cat(" Average Job Reach Using Only Formal Modes :", GCR_avg_jobs_formal, "\n")
Average Job Reach Using Only Formal Modes : 168332.3
# # Average accessibility using formal Modes (%):
GCR_avg_acc_formal <- as.numeric((cairo_hexagons$access_per_60_formal %*% cairo_hexagons$pop2018cap) / sum(cairo_hexagons$pop2018cap))
cat(" Average Accessibility Using Only Formal Modes :", GCR_avg_acc_formal)
Average Accessibility Using Only Formal Modes : 2.875381
Let’s check the regional level scores (scores for central, inner and outer Cairo)
library(dplyr)
cairo_df <- as.data.frame(cairo_hexagons) #convert sf to dataframe
cairo_df %>% group_by(znng_t_) %>%
summarise(score_formal = (access_per_60_formal %*% pop2018cap) / sum(pop2018cap), # formal transit modes
jobs_formal = (jobs_60_formal %*% pop2018cap) / sum(pop2018cap),
jobs_all = (jobs_60_all %*% pop2018cap) / sum(pop2018cap),
score_all = (access_per_60_all %*% pop2018cap) / sum(pop2018cap), # all transit modes
inc_abs = (score_all - score_formal), # accessibility increase
inc_rel = ((score_all - score_formal)/ score_formal)*100, # % increase FROM formal modes TO all modes
inc_factor = (score_all/score_formal)) # increase as a multiple
The analysis has now been run twice: 1) All public Transport Modes 2) Formal Public Transport Modes
Let’s check the difference in job reach
# Add a column showing the difference in job reach (i.e how much informal transit extends job reach)
cairo_hexagons$jobs_60_inf = cairo_hexagons$jobs_60_all - cairo_hexagons$jobs_60_formal
The above gives negative values for 12 hexagons. These values are all close to 0. This is an otp issue since it does not make sense for isochrones to shrink when more transit options are available. To avoid negative values, I use a for loop to replace them with zeros
for (i in 1:nrows){
if (cairo_hexagons$jobs_60_all[i] < cairo_hexagons$jobs_60_formal[i]){
cairo_hexagons$jobs_60_inf[i] = 0
}
else{
cairo_hexagons$jobs_60_inf[i] = cairo_hexagons$jobs_60_all[i] - cairo_hexagons$jobs_60_formal[i]
}
}
Difference in accessibility scores
# Get additional % of jobs reached.
cairo_hexagons$access_per_60_inf = cairo_hexagons$access_per_60_all - cairo_hexagons$access_per_60_formal
We also have the issue of negative values for the same 12 hexagons. For loop it is:
# The above will also give 12 negative values (same issue as jobs_60_inf). For loop used:
for (i in 1:nrows){
if (cairo_hexagons$access_per_60_all[i] < cairo_hexagons$access_per_60_formal[i]){
cairo_hexagons$access_per_60_inf[i] = 0
}
else{
cairo_hexagons$access_per_60_inf[i] = cairo_hexagons$access_per_60_all[i] - cairo_hexagons$access_per_60_formal[i]
}
}
Lets visualize how accessibility varies throughout the study area
library(sf)
# to plot metro as line
cairo_metro <- st_read("Cairo Shapefiles/Metro_Trips_TfC.shp")
Reading layer `Metro_Trips_TfC' from data source `/Users/Hussein/Documents/Code Repos/GIS-Coursework/Cairo Shapefiles/Metro_Trips_TfC.shp' using driver `ESRI Shapefile'
Simple feature collection with 6 features and 23 fields
geometry type: LINESTRING
dimension: XY
bbox: xmin: 31.19799 ymin: 29.84915 xmax: 31.34906 ymax: 30.16357
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
# to add text labels for the new towns on the outskirts
cairo_new_towns <- st_read("Cairo Shapefiles/Central-Inner_Shiyakha_NDC_CAPMAS.shp") %>%
filter(zoning_tfc == "Outer")
Reading layer `Central-Inner_Shiyakha_NDC_CAPMAS' from data source `/Users/Hussein/Documents/Code Repos/GIS-Coursework/Cairo Shapefiles/Central-Inner_Shiyakha_NDC_CAPMAS.shp' using driver `ESRI Shapefile'
Simple feature collection with 653 features and 18 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 30.86378 ymin: 29.7111 xmax: 31.83973 ymax: 30.35001
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
# Remove 10th of Ramadan City and Giza_Outer Labels)
cairo_new_towns <- cairo_new_towns[!duplicated(cairo_new_towns$name_citya),] %>%
filter(!grepl('Giza_Outer|10th of Ramdan', name_citya))
library(tmap)
tmap_mode("plot")
tmap mode set to plotting
# breaks argument used instead of style
breaks = c(0, 10000, 20000, 100000, 500000, 1000000, 2500000)
tm_shape(cairo_hexagons) +
tm_fill("jobs_60_all",
#style = "jenks", # used instead of user defined breaks
breaks = breaks,
palette = 'GnBu', # for specific colors: c('#d7191c', '#fdae61', '#ffffbf', '#abdda4', '#2b83ba', '#253494')
#legend.hist = TRUE,
title = "Jobs Within 60 min \n(Public Transit, AM Peak)") +
tm_layout(title = "Accessibility Across the GCR", # add a title
title.size = 1.5,
title.color = "azure4",
title.position = c("left", "top"),
inner.margins = c(0.09, 0.10, 0.10, 0.08), # increase map margins to make space for legend
fontfamily = 'Georgia',
#bg.color = "grey95",
frame = TRUE) +
tm_borders(col = "grey40", lwd = 0.1)+
tm_legend(title.size=0.9,
text.size = 0.6,
#frame = "grey",
position = c("right", "bottom")) +
tm_scale_bar(color.dark = "gray60",
position = c("left", "bottom")) +
tm_shape(cairo_metro) +
tm_lines(col = 'firebrick4', lwd = 1.5, alpha = 0.8) +
tm_add_legend(type = "line", labels = 'Cairo Metro', col = 'firebrick4', lwd = 1.5) +
tm_shape(cairo_new_towns) +
tm_text(text = "name_citya", col = 'black', size = 0.7,
alpha = 0.7, bg.color = "white", bg.alpha = 0.5)
tm_shape(cairo_hexagons) +
tm_fill(c("jobs_60_all", "jobs_60_formal"),
#style = "jenks", # used instead of user defined breaks
breaks = breaks,
palette = "GnBu",
#legend.hist = TRUE,
title = "Jobs Reachable \nin 60 min (AM Peak)") +
tm_layout(title = c("All Public Transit", "Formal Public Transit"), # add a title
title.size = 1.2,
title.position = c("left", "top"),
title.color = "azure4",
inner.margins = c(0.09, 0.10, 0.10, 0.08), # increase map margins to make space for legend
fontfamily = 'Georgia',
#bg.color = "antiquewhite",
frame = TRUE) +
tm_borders(col = "grey40", lwd = 0.1)+
tm_legend(title.size=0.7,
text.size = 0.5,
#frame = "grey",
position = c("right", "bottom")) +
tm_scale_bar(color.dark = "gray60",
position = c("left", "bottom")) +
tm_shape(cairo_metro) +
tm_lines(col = 'firebrick4', lwd = 1, alpha = 0.8) +
# add legend for the metro
tm_add_legend(type = "line", labels = 'Cairo Metro', col = 'firebrick4', lwd = 1) +
# two plots together on one row
tm_facets(nrow = 1)
breaks_diff = c(0, 100, 20000, 100000, 500000, 1000000)
tm_shape(cairo_hexagons) +
tm_fill("jobs_60_inf",
#style = "pretty", # used instead of user defined breaks
breaks = breaks_diff,
palette = "OrRd",
#legend.hist = TRUE,
title = "Additional Jobs Reached \nWhen Using Informal Transit") +
tm_layout(title = "Effect of Informal Transit on Accessibility", # add a title
title.size = 1.2,
title.color = "azure4",
title.position = c("left", "top"),
inner.margins = c(0.09, 0.10, 0.10, 0.08), # increase map margins to make space for legend
fontfamily = 'Georgia',
#bg.color = "grey95",
frame = FALSE) +
tm_borders(col = "grey40", lwd = 0.1)+
tm_legend(title.size=0.9,
text.size = 0.6,
#frame = "grey",
position = c("right", "bottom")) +
tm_scale_bar(color.dark = "gray60",
position = c("left", "bottom")) +
tm_shape(cairo_metro) +
tm_lines(col = 'gray23', lwd = 1.5, alpha = 0.8) +
# add legend for the metro
tm_add_legend(type = "line", labels = 'Cairo Metro', col = 'gray23', lwd = 1.5)